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Complex networks are an important paradigm of modern complex systems sciences which allows 
quantitatively assessing the structural properties of systems composed of different interacting 
entities. During the last years, intensive efforts have been spent on applying network-based con- 
cepts also for the analysis of dynamically relevant higher-order statistical properties of time 
series. Notably, many corresponding approaches are closely related with the concept of recur- 
rence in phase space. In this paper, we review recent methodological advances in time series 
analysis based on complex networks, with a special emphasis on methods founded on recurrence 
plots. The potentials and limitations of the individual methods are discussed and illustrated 
for paradigmatic examples of dynamical systems as well as for real-world time series. Complex 
network measures are shown to provide information about structural features of dynamical sys- 
tems that are complementary to those characterized by other methods of time series analysis 
and, hence, substantially enrich the knowledge gathered from other existing (linear as well as 
nonlinear) approaches. 
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1. Introduction 

The understanding of principles and mechanisms underlying the dynamics of natural systems is closely 
related to the progress of complex systems analysis. Concepts originated in the field of nonlinear dynam- 



ics such as correlation dimension |Grassberger Procaccia 1983| or Lyapunov exponents [Wolf et al. 



1985 have been introduced and successfully applied for quantitatively describing phase space topology 



and resulting dynamical properties. Spatially extended systems have been studied using, e.g., fractal prop- 



*Electronic address: reik.donner@pik-potsdam.de 
^Electronic address: ensmall@polyu.edu.hk 
^Electronic address: donges@pik-potsdam.de 
§ Electronic address: marwan@pik-potsdam.de 
^Electronic address: yong.zou@pik-potsdam.de 
"Electronic address: ms .ruoxi-xiang@polyu. edu.hk 
"Electronic address: juergen.kurths@pik-potsdam.de 



1 



2 



erties Marwan et al. , 2007c Dombradi et al. , 20071 , information-based measures Schirdewan et al. , 2007] 



or complex network approaches [Donges et al, 2009a I. 



In the past two decades, a new class of dynamical characteristics has received increasing attention, 
which is based on the widely observed phenomenon of recurrences |Marwan et al. 2007b . Many dynamical 
processes exhibit recurrences, which have already been recognized by Poincare |1890 in his seminal study 
of the three-body problem. In the context of time series analysis, we refer to a recurrence of a state Xi at 
time t = i-At (where i £ N, At is the sampling time, and x £ W 71 a state in the m-dimensional phase-space) 
whenever the state of the system Xj at another time j ■ At is similar to that initial state (i.e., Xi ~ Xj) 
or as close as we wish (but usually not identical)^] Despite the implicit technical restriction of constant 
sampling time made in this definition, we would like to note that unlike some other basic approaches of 
time series analysis, the recurrence concept can be directly generalized to unequally sampled data. 

The increasing interest in using the concept of recurrence for the analysis of dynamical systems is 



related to the introduction of more and more powerful computers Marwan, 2008 1. First return maps and 



recurrence time statistics have been introduced to study chaotic dynamical systems, unstable periodic 
orbits, or dynamical invariants |Procaccia et al. 1987; Gao, 1999 . Eckmann et al. 1987| have introduced 
recurrence plots (RPs) for visualization of recurrences in phase space. A RP represents all recurrences in 
form of a binary matrix R, where Rij = 1 if the state Xj is a neighbor of X{ in phase space, and R4 a = 
otherwise. 

RPs can be defined in different ways. In the original RP definition of Eckmann et al. |1987 , only the 
k nearest neighbors of states in phase space are considered. This preserves a constant column sum in R, 
i.e., the recurrence point density (or local recurrence rate) 



1 N 



(1) 



is conserved at RRi = k/N (with N being the length of the time series). The advantage of this method 
is that it allows comparing RPs of different systems without the necessity of normalizing the underlying 
time series beforehand, since the global recurrence rate 

1 N 

RR = JriYl Ri >i ( 2 ) 



is fixed at the same value. Alternatively, in the most common definition of a RP, a state is considered to 
be recurrent if the system's trajectory approaches state Xi in phase space closer than a certain recurrence 
threshold e, i.e., 

Rij(e) = Q(e -\\xi- Xj\\) , (3) 

where O(-) is the Heaviside function and ||-|| is a norm. The basic principle is illustrated in Fig. [l]for one 
realization of the Lorenz system 




a(y - x) 
x(r - y) 

xy — f3z 



(4) 



Further definitions of recurrences add dynamical aspects, such as local rank orders or strictly parallel 
evolution of states (parallel segments of phase-space trajectory considered in iso-directional RPs |Horai 
eTar|[2002] ). For a more detailed overview, we refer to [Marwan et al.\ |2007bl iBandt et al.\ 12008 



RPs of dynamical systems with different types of dynamics exhibit distinct structural properties (see 
Fig. [2), which can be characterized in terms of their associated small-scale as well as large-scale fea- 
tures Marwan et al. , 2007b . A periodic regime is reflected by long and non-interrupted diagonal lines. 



1 Whenever we refer to a state in phase-space, we consider either a system with all known system variables, or a phase space 
which is reconstructed from a time series, e.g., by means of time-delay embedding Packard et al. 1980 Takens 1981 . In the 
following, we presume that the reader is familiar with embedding techniques and estimation of embedding parameters. 
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Fig. 1. Basic concepts beyond recurrence plots and the resulting recurrence networks, exemplified for one realization of the 
Lorenz system (Eq. Q) with the parameters r = 28, a = 10 and /3 = 8/3 (sampling time At — 0.02, original coordinates, 
no embedding, recurrences defined based on a fixed threshold e = 5.0 using maximum norm). (A) A state at time i (red 
dot) is recurrent at another time j (black dot) when the phase space trajectory visits its close neighborhood (gray circle). 
This is marked by value 1 in the recurrence matrix at States outside of this neighborhood (small red circle) are marked 

with in the recurrence matrix. (B) Graphical representation of the corresponding recurrence matrix (recurrence plot) and 
adjacency matrix (modulo main diagonal). (C) A particular path in the recurrence network for the same system embedded in 
the corresponding phase space. 

The vertical distance between these lines corresponds to the period of the oscillation. A chaotic dynamics 
also leads to diagonals, which are however clearly shorter. There are also certain vertical structures, which 
are not as regular as in the case of a periodic motion. The RP of an uncorrelated stochastic signal consists 
of mainly isolated black points. The distribution of the points in such a RP appears rather erratic but 
nevertheless homogeneous. 
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Fig. 2. Exemplary recurrence plots of (A) a periodic motion with one frequency, (B) the chaotic Lorenz system (same 
trajectory as in Fig. |TjA) , and (C) of normally distributed white noise. 




The study of recurrences by means of RPs has become popular with the introduction of recurrence 
quantification analysis (RQA) |Zbilut &z Webber Jr. 1992 Marwan et al. , 2002b . The initial purpose of this 
framework has been to introduce measures of complexity which distinguish between different appearances 
of RPs Marwan 2008] , since they are linked to certain dynamical properties of the studied system. RQA 



measures use the distribution of small-scale features in the RP, namely individual recurrence points as well 
as diagonal and vertical line structures. RQA as a whole has been proven to constitute a very powerful tech- 
nique for quantifying differences in the dynamics of complex systems and has meanwhile found numerous 



applications, e.g., in astrophysics |Kurths et al. 


1994 , ecology Facchini et al. 


2007 , engineering 


et al. 


2009 , geo- and life sciences [Marwan et al. 


2003, 


2007a , or protein research [Giuliani et al. 



Litak 



2002 



4 



Zbilut et al. , 2004 . For a more comprehensive review on the potentials of this method, we refer to Marwan 



Zbilut et al. 


2004 . For ; 


2008 


Webber Jr. et al. 



2009 



In addition, we would like to remark that even dynamical invariants, like 
the K<i entropy and mutual information, or dimensions (information and correlation dimensions D\, D2) 
can be efficiently estimated from RPs Thiel et al., 2004; Marwan et al., 2007b . Moreover, RPs have also 



been successfully applied to study interrelations, couplings, and phase synchronization between dynamical 



2010 



. 2002a 


Romano et al. 


2004 


2005, 


2007 


Van Leeuwen et al. 


2009 


Nawrath et al. 



Another appealing concept for analyzing structural features of complex systems is based on their repre- 
sentation as complex networks of passive or active (i.e., mutually interacting) subsystems. An undirected, 
unweighted complex network G, consisting of ./V vertices and E edges, is conveniently represented by the 
binary adjacency matrix A, where Aij = 1 if vertex i connects to vertex j, and A^a = if the edge 
does not exist. 

Starting from mathematical results on graph theory, numerous applications of complex networks have 



been considered in the literature, including studies of networked infrastructures Amaral et al, 2000 La- 



social interactions Freeman, 1979 



tora & Marchiori, 2001; Guimera et al., 20051, the derivation of network patterns from empirical data of 



the assessment of functional connectivity in the brain from spatially 

or the identifica- 



distributed (multi-channel) neurophysiological measurements Zhou et al. , 2006 , 2007 



tion of dynamically relevant backbone structures in complex network representations of continuous systems 
such as atmospheric dynamics |Donges et al. 2009b|a , to mention only some important recent fields of 
application. For a more detailed statistical description of the topological features of real-world as well as 
model networks, a large variety of different statistical measures have been suggested |Albert fe Barabasi 



2002 Newman, 2003 Costa et al. , 2007 . These measures have been successfully applied to quantify the 



properties of complex networks in various scientific disciplines, fostering substantial progress in our under- 



standing of the interplay between structure and dynamics of such networks Wang Sz Chen , 2002 ; Boccaletti 



et al. , 2006 Arenas et al. 2008 



We emphasize that there are strong conceptual similarities between, on the one hand, the reconstruction 
of network topologies from spatially distributed time series (e.g., in neurophysiological or climate networks) 
and, on the other hand, the study of phase space properties of dynamical systems based on individual time 
series. Following this idea, fundamental characteristics of a dynamical system can be captured by properly 
defining complex networks based on such time series. Among other methods, the re-interpretation of the 
recurrence matrix R as the adjacency matrix A of an unweighted complex network (Fi gs. [I] and [3|) provides 
a novel concept for nonlinear time series analysis [Marwan et al. 2009 Donner et al.\ 2010b|a . 




Fig. 3. A graphical representation of the Lorenz attractor based on the recurrence matrix represented in Fig. [TJ The color of 
the vertices corresponds to their temporal order (from orange to bright green). 



The remaining part of this paper is organized as follows: In Sec. [2j we review recent approaches for 
transforming time series into complex network representations. Specifically, complex networks based on 
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Table 1. Summary of the definitions of vertices and the criteria for the existence of edges in existing complex network approaches 
to time series analysis. 



Method 


Vertex 


Jidge 


Directedness 


Proximity networks 








Cycle networks 


Cycle 


Correlation or phase space distance between cycles 


undirected 


Correlation networks 


State vector 


Correlation coefficient between state vectors 


undirected 


Recurrence networks 








k-nearest neighbor networks 


State (vector) 


Recurrence of states (fixed neighborhood mass) 


directed 


adaptive nearest neighbor networks 


State (vector) 


Recurrence of states (fixed number of edges) 


undirected 


e-recurrence networks 


State (vector) 


Recurrence of states (fixed neighborhood volume) 


undirected 


Visibility graphs 


Scalar state 


Mutual visibility of states 


undirected 


Transition networks 


Discrete state 


Transitions between states 


directed 



different definitions of RPs (so-called recurrence networks) are discussed in some detail. Section [3] summa- 
rizes technical issues that have to be considered when systematically applying the different approaches to 
time series analysis. Finally, two examples for real- world applications of recurrence networks are discussed 
in Sec. HI 



2. Transforming time series into complex networks 

Recently, several approaches have been proposed for transforming (observational) time series into complex 
network representations. These methods can be roughly distinguished into three classes (see Tab. [I]), which 
are based on 

(i) mutual proximity of different segments of a time series (proximity networks), 

(ii) convexity of successive observations (visibility graphs), and 

(hi) transition probabilities between discrete states (transition networks). 

With the exception of visibility graphs, all approaches are related with the concept of recurrence. This 
is particularly evident for proximity networks, where connectivity is defined in a data- adaptive local way, 
i.e., by considering distinct regions with a varying center at a given vertex in either the phase space itself 
or an abstract proximity space. In contrast, for transition networks, the corresponding classes are rigid, 
i.e., determined by a fixed coarse-graining of the phase space. The distinction between both classes of 
approaches is conceptually similar to the duality of symbolic time series analysis (i.e., time series analysis 



based on a coarse-graining of the dynamics) and quantitative analysis of RPs |Donner et al. 2008 , which 



may both be used for estimating similar dynamical invariants such as entropies and mutual information. 

In addition to these specific relationships between the recurrence concept and different types of time 
series networks, there is a fundamental structural analogy between RPs and (unweighted) complex net- 
works in general. Both structures are based on binary matrices (i.e., recurrence and adjacency matrices, 
respectively) that can be used for studying basic topological properties of the underlying complex system 
based on sophisticated statistical measures. Proximity and transition networks as well as RPs based on 
Eq. (§ can be generalized by withdrawing the application of a specific threshold, which leads to weighted 
networks and unthresholded RPs (distance plots), respectively. For example, the unthresholded RP ob- 
tained from one trajectory of a given dynamical system may be re-interpreted as the connectivity matrix 
of a fully coupled, weighted network. 

Among the three classes of methods listed above, the largest group of concepts is given by proximity 
networks, where the mutual closeness or similarity of different segments of a trajectory can be characterized 
in different ways. Consequently, there are different types of such proximity networks (see Tab. [TJ: cycle 
networks, correlation networks, and recurrence networks. However, all these methods are characterized by 
two common general properties: 

Firstly, the resulting networks are invariant under relabeling of their vertices in the adjacency matrix. 
Hence, the topological characteristics of proximity networks yield nonlinear measures that are invariant 
against permutation of vertices. In this respect, the network-theoretic approach is distinctively different 
from traditional methods of time series analysis where the temporal order of observations does explicitly 
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matter. 

Secondly, we have to point out that particularly proximity networks are spatial networks. In particular, 
recurrence networks are embedded in the phase space of the considered system, with distances being defined 
by one of the standard metrics (e.g., Euclidean, Manhattan, etc.). Similar considerations apply to other 
types of proximity networks as well. 

Both mentioned characteristics imply that the network-theoretic concept of a path on a given graph (see 
Fig. [T]C) is distinctively different from the trajectory concept that records the causal dynamic evolution 
of the system [Donner et al. 2010b] . Note that unlike for proximity networks, causal relationships are 



conserved in transition networks (and at least to some extent also in visibility graphs). 

In the following, we will discuss the main properties of the different concepts in some detail. 



2.1. Cycle networks 

Zhang & Small |2006j (see also |Zhang et al. 2008, Small et al. 2009 1) first suggested to study the topological 



features of pseudo-periodic time series by means of complex networks. Suppose that a dynamical system 
possesses pronounced oscillations (examples are the well-known Lorenz and Flossier systems). In this case, 
we identify the individual cycles contained in a time series of this system with the vertices of an undirected 
network. Edges between pairs of vertices are established if the corresponding segments of the trajectory 
behave very similarly. For quantifying the proximity of cycles in phase space, different measures have been 



proposed. Zhang et al. 2006 introduced a generalization of the correlation coefficient applicable to cycles 



of possibly different lengths. Specifically, this correlation index is defined as the maximum of the cross 
correlation between the two signals when the shorter of both is slid relative to the longer one. That is, if 
the two cycles being compared are C\ = {x\, x 2 , ■ ■ ■ , x a } and C2 = {2/1, 2/2, • • • , Up] with (without loss of 
generality) a < (3, then we compute 



p(Ci,C 2 )= max {(x 1 ,x 2 ,...,x a ),(yi + i,y2+i,..-,y a +i)) , 

i=a,...(/3-a) 

where (•, •) denotes the standard correlation coefficient of two a-dimensional vectors, and set 



A 



1, .1 



@(p(Ci, Cj) - Pmax) ~ 5, 



':)■ 



(5) 



(6) 



where <5jj is the Kronecker delta necessary in order to obtain a network without self-loops. As an alternative, 
the phase space distance Zhang & Small, 2006 



D(d,C 2 ) 



mm 



i=0,...((3-a) a 



1 a 



.X j 



Vj+il 



has been suggested, leading to the following definition: 



A 



1 j 



D(Ci,Cj))-5, 



(7) 



(8) 



Of course, there are other calculations one could perform as well. 

As an example for constructing complex networks from time series, Fig. [4] shows one realization of the 
Lorenz system, which is characterized by a double-scroll topology of the attractor with pronounced chaotic 
oscillations. For the first about 10 time units of simulation, the system rotates around one of both unstable 
centers (x, y < 0), is then captured by the other part of the attractor (x, y > 0) for about 5 more time units, 
followed by some fast transitions between both parts (i.e., involving only one or two subsequent oscillations 
around each center), before performing again rotations around the second center (x, y > 0) between t ~ 25 
and t ~ 32. This structure is well represented in the adjacency matrix A of the corresponding cycle network 
based on the x-coordinate time series (see Fig.[5j^), where we observe several pronounced clusters along the 
main diagonal corresponding to the two distinct parts of the attractor. Consequently, the resulting network 
(Fig. [6]A) shows a pronounced community structure with two groups, corresponding to the double-scroll 
topology of the system. 

The advantage of cycle networks is that explicit time delay embedding is avoided. In addition, the 
method is more robust than other methods against additive noise, given a small enough noise magnitude 
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Fig. 4. Time series of the trajectory of the Lorenz system (Eq. Q) used in Figs. [T] [2] and [3] (At = 0.05). 



to allow a clear identification of the individual cycles from the time series. Moreover, cycle networks are 
invariant under reordering of the cycles (this is precisely the same property that was also exploited for cycle- 
shuffled surrogate methods [Theiler fc Rapp 1996| but not the pseudo-periodic surrogate method Small 
et al., 2001] ) . However, for chaotic and nonlinear systems in a near-periodic regime, we typically observe 
significant orderly variation in the appearance of individual cycles. For systems that are linear or noise 
driven, that orderly variation will be less pronounced. As a consequence, the networks constructed with 
these methods will have characteristic and distinct properties: linear and periodic systems have cycle 
networks that appear randomly, while chaotic and nonlinear systems generate highly structured networks 
Zhang & Small, 2006 Zhang et al. , 2008 . Therefore, the vertex and edge properties of the resultant 



networks can be used to distinguish between distinct classes of dynamical systems. Moreover, Zhang & 
Small 2006 used meso-scale properties of the networks — and in particular the clustering of vertices — to 
locate unstable periodic orbits (UPOs) within the system. This approach is feasible, since a chaotic system 
will exhibit a dense hierarchy of unstable periodic orbits, and these orbits act as accumulation points in 
the Poincare section. Hence, the corresponding vertices form clusters in the cycle network. 



2.2. Correlation networks 



By embedding an arbitrary time series, individual state vectors Xi in the m- dimensional phase space of 
the embedded variables can be considered as vertices of an undirected complex network. Specifically, if the 
Pearson correlation coefficient 



(9) 



is larger than a given threshold r, the vertices i and j are considered to be connected Yang & Yang, 2008 



Gao & Jin 2009a 



Aij = 6(r 



(10) 



Interpreting 1 



as a proximity measure, the condition r j ~ > r corresponds to the definition (13 



of a recurrence with e = 1 



r. 



The consideration of correlation coefficients between two phase space 



vectors usually requires a sufficiently large embedding dimension m for a proper estimation of rij. Hence, 
information about the short-term dynamics might get lost. Moreover, since embedding is known to induce 



spurious correlations Thiel et al., 2006 , the results of the correlation method of network construction may 



suffer from related effects. 

The adjacency matrix of a correlation network is shown in Fig. [5^3 for the same example trajectory of 
the Lorenz system as previously used for constructing a cycle network. In contrast to the cycle network, we 
observe strong connectivity between vertices corresponding to the time intervals t = . . . 7 and t = 24 ... 28, 
i.e., two time intervals where the trajectory is actually captured in two different parts of the attractor. An 
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Fig. 5. Adjacency matrices corresponding to different types of networks constructed from the ^-coordinate of the Lorenz 
system shown in Fig. [4] (A) Cycle network (JV = 40, critical cycle distance in phase space D max — 5), (B) correlation network 
(TV — 654, embedding dimension m — 10 with delay r = 3 time steps), (C) fc-nearest neighbor network (asymmetric version), 
TV = 675, m = 3, r = 3, k = 10, corresponding to a recurrence rate of RR ~ 0.015 using Euclidean norm; the associated 
adaptive nearest neighbor network (not shown) is characterized by a very similar pattern), (D) £-recurrence network (JV = 675, 
m = 3, r = 3, e = 2, maximum norm), (E) visibility graph (TV = 681), and (F) transition network (based on an equipartition 
of the range of observed values into TV = 20 classes of size Ax = 3.0, minimum transition probability p = 0.2 during 3 time 
steps). Note that only in panels (C) and (D), the adjacency matrices correspond to recurrence matrices of the underlying 
time series according to the standard definition Eckmann et al. 1987 Marwan et al. 2007b . In both cases, recurrence points 



originated from strong tangential motion (sojourn points) have been removed, resulting in additional asymmetries. 



explanation for this behavior is that the dynamics itself within the two time intervals appears to be rather 
similar (see Fig. [4]), but it is just shifted in the x- (and y-) coordinate. Since the estimation of correlation 
coefficients between embedding vectors explicitly removes the mean position of the trajectory during the 
different time intervals covered by these vectors, the two respective parts of the trajectory are considered 
to be similar with respect to the correlation criterion, although they are actually well separated in the 
actual phase space. This example underlines that correlations must be carefully distinguished from true 
metric distances. 



<) 




Fig. 6. Graphical representation of the different complex networks based on the adjacency matrices shown in Fig. [5| The 



graphs have been embedded into an abstract two-dimensional space using a force directed placement algorithm Battista et al. 



1994 . For panels (A)-(E), the vertex color indicates the temporal order of observations (from orange to bright green), for the 
transition network (panel (F)), colors correspond to the different x values. Note that in panels (B) and (D), some individual 
disconnected vertices have been removed from these network representations. 



Visualization of the correlation network embedded in an abstract two-dimensional space (Fig. [6j3) 
reveals a pronounced community structure with two major groups that are characterized by a ring-like 
topology. However, these two groups do not correspond to the two scrolls of the attractor, as is the case for 
the cycle network and also for recurrence networks (see Sec. 2.3). In contrast, it is a reasonable assumption 
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that the observed group structure is determined by the orientation of the arc-like embedding vectors around 
each of the centers. 

2.3. Recurrence networks 

A recurrence network is a complex network whose adjacency matrix is given by the recurrence matrix of a 
time series^J i.e., we define the adjacency matrix of a recurrence network by 

Aij = Rij — 5ij. (11) 

Note that removing the line of identity from the RP corresponds to the consideration of the smallest 



possible Theiler window in traditional RQA | Marwan et al. , 2007b 



Since information about the temporal ordering of observations is not explicitly regarded in a recurrence 



network defined according to Eq. (11), the topological features of the resulting graphs reflect dynamically 
invariant properties associated with the specific dynamical system. From this perspective, the quantitative 
analysis of recurrence networks, although being based on the same recurrence matrix as traditional RQA, 
reflects distinctively different properties of the system than line-based RQA measures. Hence, besides RQA 
and the estimation of dynamical invariants based on line structures in RPs, the analysis of recurrence 
networks can be considered as a third column for the quantitative recurrence-based characterization of 
phase space properties of dynamical systems. Moreover, while the appropriate estimation of most RQA 
measures requires the careful choice of a second parameter (the minimum line length l m i n ), quantitative 
characteristics of recurrence networks involve only a single parameter (depending on the specific algorithm, 



see below). However, computing network-theoretic measures (e.g., betweenness centrality) |Newman 2003 
often requires larger computational efforts than traditional RQA. 

Since the recurrence matrix can be defined in different ways (see Sec. [I]), there are distinct sub- types 
of recurrence networks that are characterized by somewhat different structural properties: 

2.3.1. k-nearest neighbor networks 

Following the original definition of a RP by Eckmann et al. 1987|, every (possibly embedded) observation 



vector is considered as a vertex i, which is then linked to those k other vertices j that have the shortest 
mutual distances dij with respect to i in phase space (i.e., to its k nearest neighbors). This means that a 

r(k) r(k) 

directed edge is introduced from % to every vertex j G , where A/^ is the set of k nearest neighbors of 
i (see Tab. [2J. Hence, the neighborhoods defined in this way preserve a constant mass (i.e., the number of 
vertices is the same in all neighborhoods). Unlike for cycle and correlation networks, the adjacency matrix 
of the k-nearest neighbor network defined in such a way is generally asymmetric, since j £ A^f does not 

r(k) 

imply i S JVj ■ Hence, the resulting networks are characterized by directed edges. Note that an undirected 



and symmetric version can easily be obtained by setting Rj^ = 1 whenever Rij = 1 |Shimada et al. 



|2008| . 

For a fc-nearest neighbor network, the distribution of out-degrees is always fixed at P(k out ) = 5(k). In 
contrast, the distribution of in-degrees allows for some variability, but must necessarily have a mean value 
'h™} = k, since there are exactly Nk directed edges by definition (note that transforming the k- nearest 



neighbor network into an undirected graph [Shimada et al. 2008] leads to a network with Nk/2 to Nk 



undirected edges). The remaining spatial pattern of in-degrees provides information about the attractor 
geometry. Specifically, we can infer that if k l ™ <C k, v lies in a phase space region with decreased density 
compared to the surrounding attractor. In contrast, if S> k, v must be located in a densely populated 
region of the attractor. 

If the coordinates of the individual vertices in the underlying phase space are known, we can consider 
the neighborhood size 

l™ ax (k) = max^j-llxi - Xj\\} = max {\\xi - xj\\} (12) 



2 In 



Donner et al. 



2010b 



the term recurrence network has been more specifically used for the special case in which 



recurrences have been defined according to p|) 
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Table 2. Comparison of algorithms used for constructing k- and adaptive nearest neighbor networks, where 
steps (i)-(iii) are identical for both algorithms. Si is an ordered set initially containing all the nearest neighbor 
indices of i in increasing order of phase space distance and excluding i itself, i. e. , Vn < w : D i g. ^ < D i sdw)- 
Si\v denotes the removal of nearest neighbor index v from the set Si, hence, Si(l) gives the index of the 
closest neighbor remaining in the set. For the adaptive nearest neighbor networks, a simplified algorithm is 
presented where the vertices are processed in temporal order, i.e., starting from the earliest and ending with 
the latest time index, as was the convention in |Xu et al\ |2008] . The algorithm is readily generalized to an 
arbitrary processing order. 

fc-nearest neighbor network Adaptive nearest neighbor network 

(i) Calculate distance matrix D(i, j). 

(ii) Obtain Si by implicitly sorting ith row of D(i,j). 

(iii) Initialize adjacency matrix: A(i,j) = V (i,j). 

(iv) Fill adjacency matrix: 

FOR i G {1, . . . ,N} FOR j G {1, . . . ,E } 

FOR j G {1, . . . , k} FOR i G {1, . . . , N} 

v = Si(l) v = Si(l) 

A(i,v) = l A(i,v) = A(v,i) = 1 

Si = Si \ v Si = Si \ v 

Sy = Sy \ j 

as a measure that is directly related with the inverse state density p(xi)~ l of the system in the vicinity 
of a vertex i. Prom a statistical perspective, this strategy for retaining information about the attractor 
geometry can be regarded as a kernel density estimate with a simple constant kernel function, where k 
serves as a smoothing parameter (small k: good spatial resolution, but large variance of the estimated 
state density; large k: small variance, but bad spatial resolution). Note that the degree centrality of an 



e-recurrence network (see Sec. 2.3.3) can be interpreted in a similar way, with the neighborhood size e 
serving as the smoothing parameter. 

Figure [5p displays the adjacency matrix of a k- nearest neighbor network for the Lorenz system, which 
corresponds to the respective RP (modulo the main diagonal). We also note the strong similarity with 
the connectivity of the associated cycle network. Reembedding the network graphs into a two-dimensional 
space (Fig. [6p) allows recovering the double-scroll pattern of the original attractor in the reconstructed 
phase space of the three-dimensional embedding vectors. Note that the community structure with two 
ring-like network components actually reflects the different parts of the attractor. 

2.3.2. Adaptive nearest neighbor networks 

Unlike other approaches for transforming time series into complex networks, the fe-nearest neighbor method 
leads to directed networks. However, in many cases the properties of undirected networks would be more 
directly interpretable. Moreover, the total number of undirected edges E is not fixed by the algorithm 
itself. Specifically, there are some vertices with k™ < k, which has certain disadvantages if one wishes to 
study, e.g., the distributions of motifs (i.e., small subgraphs consisting of a fixed, low number of vertices) 
of a given order in the network. 



In order to define an undirected nearest neighbor network with a precise control of E, Xu et al. 



2008 1 as well as Small et al. |2009 suggested an alternative network construction method considering 



nearest neighbors but correcting for a constant number of distinct edges Eq assigned to each vertex. 
In their approach, the network construction is an iterative process (see Tab. [2]), where in each step all 
observations (vertices) are linked to their nearest neighbors in phase space. However, if vertex i is linked with 
vertex j, vertex i is removed from the neighborhood of j. This avoids the possibility of "double-counting" 
vertex i as a neighbor of vertex j and vice versa. Hence, the link between j and i is bi-directional, i.e., 
j) = A(j, i) = 1, resulting in a symmetric adjacency matrix A, (i.e., an undirected network). This edge 
construction is repeated Eq times. Finally, from each phase space vector exactly Eq edges have been drawn 
to its geometric neighbors, which thus become also neighbors in a complex network sense. Consequently, 
there are exactly NEq undirected edges, which connect vertices of at least degree Eq. Specifically, a phase 
space vector can be a neighbor of more than Eq other phase space vectors, with an average degree (k) = 2Eq. 
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In the following, we will refer to the resulting networks as adaptive nearest neighbor networks. 

The construction of adaptive nearest neighbor networks differs from the fc-nearest neighbor network, 
since the resulting matrix is symmetric, i.e., the edges defined here are undirected from the beginning. 
Nonetheless, the process is more subtle than simply symmetrizing the recurrence matrix R by taking the 
logical matrix (R + R T ) > 0. The iterative network construction method generates a matrix such that 
the closest Eq neighbors are always included. The exclusion process described above works to include the 
next closest neighbors from among the possible candidates. Note that for Eq = k, the adaptive neighbor 
network always includes all edges of the associated ^-nearest neighbor network. However, adaptive nearest 
neighbor networks always have higher edge densities than fc-nearest neighbor networks. 

The frequency distribution of motifs in adaptive nearest neighbor networks has been demonstrated 



to be a sensitive indicator of the specific type of dynamics in the underlying dynamical system. In Xu 



et al. 2008 



networks generated from various dynamical systems were compared and it was found that 
the specific distributions of the motif frequency differed qualitatively, but did so consistently. That is, 
periodic systems exhibit one particular type of distribution, but chaotic (one positive Lyapunov exponent) 
and hyper-chaotic dynamics (more than one positive Lyapunov exponent) different ones. Specifically, motif 
prevalence is determined by the heterogeneity of the attractor and the intrinsic dimensionality of the system, 
both being larger for a chaotic system than for periodic dynamics. As a consequence, non-transitive motif 
patterns are more common in chaotic systems than in periodic ones |Xu et al. , 2008] . Recently, it has been 
shown that there is a distinct relationship between the motif distributions obtained for certain stochastic 
processes, and the associated scaling exponents |Liu & Zhou 2009 . Note that other types of recurrence 



networks can be expected to show different distributions of motif frequencies. In particular, whether or not 
the motif distributions of other types of recurrence networks can be used for distinguishing qualitatively 
different dynamics as well remains a subject of future studies. 



2.3.3. e-recurrence networks 

For adaptive as well as /c-nearest neighbor networks, choosing an equal number of neighbors for each point 
in phase space allows obtaining a representation of the underlying attractor that is independent of the local 
metric properties of the attractor in the considered embedding space. Hence, the resultant networks are 
based on the relative proximity between points on a trajectory in phase space and are therefore independent 
of any monotonic rescaling of the data (as is typically permitted by the various embedding theorems in the 
guise of an observational function) — the same network will result independent of observation function. 
In theory this should be particularly useful for measuring metric invariants such as correlation dimension. 
However, it is unclear whether this method, or other types of recurrence networks are to be preferred in 
practice — presumably this distinction will depend on the particular application. 

As a disadvantage of both types of nearest neighbor networks, there is no direct relationship between 
their local as well as global properties and the invariant density of the system under study. As an alternative, 
the neighborhood of a single point in phase space can also be defined by a fixed phase space distance e 
(see Eq. [§) [Wu et al.\ |2008] |Gao fc Jin] |2009a|b] |Marwan et~al\ |2009] |Donner et al\ |2010b| , i.e., by 



considering fixed phase space volumes instead of a fixed (local or global) number of edges. In the following 
we will refer to this type of network as an e-recurrence network. Note that the structural properties of such 
networks often show a high degree of similarity with those of nearest-neighbor networks with similar link 
density (cf. Figs.[5]C!,D and[6]C,D). However, the local network properties can be directly related with the 
phase space properties of the underlying system (for a more detailed review, see |Donner et al. , 2010b]). As 
an example, Fig. [7] shows the trinity of centrality measures (degree, closeness, and betweenness |Freeman 



1979 1) as well as the local clustering coefficient for realizations of the Lorenz system at various values of 
the control parameter r: 

(i) The degree centrality k v (Fig. [7]A) gives the number of neighbors of a vertex v and is therefore pro- 
portional to the local recurrence rate RR V and, hence, the phase space density at the corresponding 
point in phase space. 

(ii) The closeness centrality c v (Fig. [7p) is related to the inverse mean network distance of a vertex with 
respect to all other vertices, implying that high values of closeness can be expected in the central parts 
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Fig. 7. (A) Degree, (B) local clustering coefficient, (C) betweenness (in logarithmic units), and (D) closeness centrality for 
the e-recurrence networks obtained from trajectories of the Lorenz system (Eq. Q) for different control parameters r. To 
enhance the visibility of the underlying structures, the networks have not directly been derived from the time series, but from 
TV = 6, 000 points (for each value of r) in the associated Poincare sections at x — 0, x < (maximum norm, e = 0.05<r with 
a being the empirical standard deviation of the considered data), z-coordinates have been suppressed in the figure. Note that 
in contrast to the standard definition of closeness Donner et at. 2010b , which is useful if only few isolated vertices exist, c v 



has been computed here separately for the individual mutually disconnected subgraphs. The broad windows in r with sparse 
points indicate the presence of periodic orbits in the system, which have not been perfectly sampled in the Poincare sections. 



of the attractor, whereas the outer parts are characterized by small values, 
(iii) The betweenness centrality b v (Fig. [Tt!) measures the number of shortest paths between pairs of vertices 
in the network that traverse a given vertex v and, hence, indicates regions of phase space that are 



characterized by low density, but separate regions of higher density (geometric bottlenecks |Donner 
et al. 2010a|). Specifically, high betweenness values indicate a strong local attractor fragmentation. 



Note that the betweenness is partially influenced be the degree: Phase space regions with high degree 
often show low betweenness, since there are many redundant shortest paths traversing this region. 
Nonetheless, betweenness centrality still yields additional information, since it is not defined exclusively 
locally, but encodes global network properties |Donner et al. 2010b|a . In particular, vertices with low 



degree, but high betweenness are of potential interest. 
We would like to remark that for a fixed e, all three centrality measures are extensive network properties 
by definition (i.e., their values depend on the system size N, either in a linear (k v ) or a nonlinear (c v , b v ) 
way). In contrast to this, the local recurrence rate RR V = k v /(N — 1) (i.e., the density of connections 
in the vicinity of a vertex v) is a non-extensive property (i.e., RR V does not depend on TV apart from 
possible finite-size effects). 

(iv) Another non-extensive vertex property is the local clustering coefficient C v (Fig. [7^>), which measures 
the presence of closed triangles in the network and, hence, characterizes localized higher-order spatial 
correlations between observations. Specifically, since recurrence networks are spatial networks, it is 
possible to interpret the structures resolved by spatial variations of C v in terms of the heterogeneity 
of the spatial filling of points. Donner et al. [2010b|a| have demonstrated that this interpretation is 



consistent with the fact that high values of C v often coincide with dynamically invariant objects, such 
as unstable periodic orbits or, more generally, invariant manifolds. For the Lorenz system, regions of 
the attractor with high C v coincide with supertrack functions in the Poincare map (corresponding the 
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unstable periodic orbits of the full system), i.e., regimes of intermittent dynamics. 



Global network measures, such as the global clustering coefficient C (i.e., 
over all vertices) , the closely related transitivity T Boccaletti et al. , 2006 
C (i.e., the mean graph distance between all pairs of vertices), are well suited for tracing qualitative 



the average value of C v taken 
and the average path length 



changes in the dynamics (see Sec. 4.2). For the global clustering coefficient, this is a direct consequence 



of the different local divergence between neighboring trajectories in case of periodic and chaotic systems. 
However, in general, we have to carefully distinguish between discrete maps and continuous dynamical 



systems Zou et al, subm. : since the topological properties of periodic trajectories in phase space strongly 



differ between both types of systems, the behavior is distinctively different between maps (small C and 
large C and T for periodic orbits, large C and small C and T for chaotic trajectories,) and continuous 
systems (opposite behavior of C for comparable RR). From this perspective, we should note that (as for 
traditional RQA), different network measures may point to similar dynamical properties (i.e., may not be 
fully independent of each other). 

A detailed discussion of the geometric interpretation of a variety of global network properties as well 
as vertex and edge properties of e-recurrence networks, including graphical representations of the spatial 
distributions of different vertex properties for the Lorenz system in the standard parameter setting (cf. 
Fig. [I]), can be found in [Donner et al. 2010b . 



2.3.4. Recurrence network analysis and RQA 

Following the above considerations, it is evident that network-theoretic measures obtained from recurrence 
networks characterize (in most cases) distinctively different properties of a complex system than RQA 
measures. Specifically, RQA measures are based on continuous line structures in recurrence plots, i.e., rely 
on temporal interdependences between individual observations (or parts of a trajectory). In contrast to this, 
temporal information is not considered in the network analysis, which therefore covers mainly geometric 
properties of the system in phase space (i.e., spatial dependences). In this respect, the nonlinear statistical 
concept that has possibly the closest similarity with recurrence network analysis is the estimation of fractal 
dimensions. However, this method is much more restrictive than the network view, since it explicitly 
assumes the presence of geometric self-similarity in phase space. Moreover, network characteristics such as 
motif distributions or clustering coefficients are based on higher-order statistical dependences, i.e., mutual 
neighborhood relationships between more than two different points in phase space. In a similar way, we 
can argue for path-based network measures (e.g., average path length or betweenness centrality). 

From these fundamental conceptual differences it follows that network-theoretic measures do indeed 
capture complementary aspects of a complex system in comparison not only to RQA, but also most other 
established methods of time series analysis. It should be noted, however, that there are certain fields of 
applications that can in principle be addressed using both traditional RQA and recurrence network analysis. 



One important example is the detection of dynamical transitions in time series (see Sec. 4.2). However, 
since both concepts provide complementary points of view, their combined use is often desirable in order 
to obtain additional information. 



2.4. Other approaches 

2.4.1. Visibility graphs 

The concept of visibility graphs as networks of intervisible locations in physical space has been known 
for decades and has found many practical usages in, among other fields, engineering and urban plan- 
ning |de Floriani et al., 1994 Turner et al. 2001 . Recently, Lacasa et al. iLacasa et al. 2008 transferred 



this concept to the field of time series analysis. Here, individual observations in a time series are identified 
with vertices of an undirected complex network, and their connectivity is established according to a local 
convexity constraint between successive observations which corresponds to a visibility condition in physical 



space. The visibility graph approach has already found various applications |Ni et al. 
20091 ILiu efal] [2UT0] [Tang L Liu[ [20091 jYang eToI 



2009; Eisner et al. 2009 Luque et al. 2009; Qian 



2009 Lacasa et al. 
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et al. , in press and is particularly interesting for certain stochastic processes where the statistical proper- 
ties of the resulting network can be directly related with the fractal properties of the time series. However 
beside the relationship between the degree distribution P(k) and the Hurst parameter of the underlying 
stochastic process Ni et al. . 2009 Lacasa et al. , 2009 , convincing links between further network-theoretic 



measures and distinct phase space properties have not been found so far. Moreover, in its currently applied 
form, the visibility graph method is restricted to univariate time series analysis. In this sense, studying the 
degree distribution of visibility graphs does not provide additional information, however, it may still have 
benefits with respect to the numerical procedures. 

As it can be seen from Fig. [5^, a visibility graph typically has a distinct topology that is characterized 
by hubs corresponding to local maxima of the considered time series. The presence of these hubs gives rise 
to a pronounced community structure, where the different network clusters reflect the temporal order of 
observations (see Fig. [6b). The mentioned general features lead to degree distributions of visibility graphs 
that are often found to be scale-free, which reflects the fractal properties of the underlying time series. 



2.4.2. Transition networks 

Coarse- graining the range of values in a time series into a suitable set of classes {S\, . . . , Sk} allows 
considering the transition probabilities tt^r = P(xi + i G Sp\xi G S„) between these classes in terms of a 
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approach is equivalent to applying a symbolic discretization with static grouping |Daw et al. , 2003 ; Donner 



et al. , 2008] to the phase space of the studied system. Unlike proximity networks, the resulting transition 
networks explicitly make use of the temporal order of observations, i.e., their connectivity represents 
causality relationships contained in the dynamics of the observed dynamical system. By introducing a 
cutoff p < 1 to the transition probability n a> p between pairs of discrete "states" S a and Sp, we obtain an 
unweighted network representation, which is, however, still directed. Note that for a trajectory that does 
not leave a finite volume in phase space, there is only a finite number of discrete "states" Si with a given 
minimum size in phase space. This implies the presence of absorbing or recurrent states in the resulting 
transition network. 

The transition probability approach is well suited for identifying such "states" (i.e., regions in phase 
space) that have a special importance for the causal evolution of the studied system in terms of betweenness 
centrality b v and related measures. However, its main disadvantage is a significant loss of information on 
small amplitude variations. Moreover, the resulting networks do not only depend on a single parameter, 
but on the specific definition of the full set of classes. Note, however, that coarse-graining might be a 
valid approach in case of noisy real- world time series, where extraction of dynamically relevant information 



hidden by noise can be supported by grouping the data |Daw et al., 2003 



In contrast to the other approaches for constructing complex networks from time series, the topology 
of transition networks depends rather sensitively on the specific choice of discretization. For the example of 
the Lorenz system (x-coordinate) shown in Figs. [5b and[6b, the resulting network pattern, however, reveals 
the spatial structure of the attractor, which is caused by the fact that for the considered relatively dense 
sampling of the trajectory, the transitions between subsequent observations in time always link regions of 
phase space that are closely neighbored. This results in the pronounced alignment of connections along 
the main diagonal of the adjacency matrix. However, within the two scrolls, there may also be transitions 
bridging several "cells" of the coarse-grained phase space, which is reflected by a higher connectivity of 
the network among the corresponding vertices. As a result, the transition network topology (represented 
in Fig. [6b as a directed graph including self-loops) again reveals the fundamental spatial structure of the 
attractor. 
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3. Practical considerations 



3.1. Cycle networks 

For an implementation of the cycle network approach, the time series must be divided into distinct cycles. In 



Zhang & Small, 2006 Zhang et al. , 2008 the preferred method for defining cycles is splitting the trajectory 



at peaks (or equally troughs). In order to quantify the mutual proximity of different cycles, different 
measures can be applied depending on the specific application. On the one hand, the cycle correlation index 
Pi j (Eq. ([5])) can be properly estimated without additional phase space reconstruction (embedding), which 
has advantages when analyzing noisy and non-stationary time series, e.g., experimental data [Zhang &| 



Small, 2006 . Moreover, this choice effectively smoothes the effect of an additive independent and identically 
distributed noise source |Zhang et al. 2006 . On the other hand, the phase space distance Dij (Eq. ([7] 



is physically more meaningful |Zhang et al. , 2008 . For the example systems as well as some real- world 



clinical electrocardiogram recordings studied in Zhang Sz Small 2006 Zhang et al. , 2008 , both methods 



have been found to perform reasonably well. However, whether the previously considered approaches also 
lead to feasible results for other cases has to be further investigated in future research. 

In general, the construction and quantitative analysis of cycle networks requires a sufficiently high 
sampling rate, i.e., we require that both cycle lengths a and j3 in Eqs. §5§ and ([7]) are reasonably large. 
The main reason for this requirement is that even two cycles that are fully identical but sampled in a 
different way may have rather different cycle correlation indices (and phase space distances) depending 
on the exact values of the observed quantity. Hence, for a very coarse sampling, it is possible that two 
cycles that are actually close in phase space may not be connected in the cycle network. However, for large 
sampling rates, the variance of this measure decreases, resulting in a more reliable network reconstruction. 



3.2. Recurrence networks 

A common problem in the construction of recurrence networks is the presence of sojourn points, which 
correspond to temporally subsequent observations within a small part of the phase space in the presence 



of strong tangential motion |Gao, 1999 . In order to avoid artificial results due to such points with strong 



temporal correlations, points belonging to the same "strand" must not be linked. As a possible solution, 



Xu et al. 1 2008 j suggested that eligible neighbors should have a temporal separation greater than the mean 



period of the data (a considerable alternative applicable also to non-oscillatory data would be the associated 
correlation time). However, removing all recurrence points with a short temporal distance can lead to a 
loss of "true" recurrences as well. Moreover, we note that this approach introduces an additional parameter 
(the minimum recurrence time). For e-recurrence networks, sojourn points can be directly removed from 



the complete recurrence matrix [Marwan et al., 2007b . For this purpose, for every vertex i, those edges 



(Aij = 1) with j > i for which j is a subsequent recurrence point of i are removed (Aij = 0). In a second 
step, the adjacency matrix is symmetrized again by setting Aji = whenever Ai j = 0. 

In any case, recurrence network properties depend on the sampling and possibly also the length of the 
time series: 

(i) For nearest neighbor networks, using a longer time series effectively leads to a finer coverage of the 
available state space. As a consequence, when keeping Eq (or k, respectively) fixed, we obtain a higher 
spatial resolution of the structural properties of nearest neighbor networks as the length N of the 
time series increases. Moreover, we note that the neighborhood of a vertex will typically change with 
increasing N since additional vertices with smaller distances in phase space appear. 

(ii) For e-recurrence networks, the explicit choice of a threshold allows to directly control the spatial 
resolution. 

In general, it is recommended to use a sampling that allows a reasonable spatial resolution of the 
whole phase space covered by the attractor. Moreover, since the choice of the reconstruction (embedding) 



parameters also matters — as for RPs Marwan et al. 2007b| — properly selected embedding parameters 
should be used wherever possible. 
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3.2.1. Adaptive nearest neighbor networks 

In contrast to /c-nearest neighbor networks, the sample network obtained from the adaptive nearest neighbor 
method will depend (slightly) on the order in which one processes the individual embedded time series points 
(originally, this was done in temporal order). Although this dependence could be completely eliminated by 
insisting that nodes are considered strictly in order of proximity to their neighbours there is no need to 
incur this additional computational complexity as the variation in the resulting network is not important. 
While there are many pathological situation under which small variation in the resultant networks can 
arise, these do not contribute to any significant structural variation in the network topology for moderate 
to large N (see Fig. [8]) . 
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Fig. 8. (A,B) Probability distribution of the Hamming distance H a p — 1 — (^S (^j 01 ) — A^) 

between the adjacency matrices of two complex networks with the same set of vertices) between different adaptive nearest 
neighbor networks obtained from the same realization of (A) the Rossler system (x = —y — z, y — x + ay, z — b + z(x — c)) 
with a — 0.15, b = 0.2 and c = 10, and (B) the Lorenz system (Eq. Q, parameters as before). In both cases, 1,000 random 
permutations of the same reference trajectory with TV = 5,000 data points have been used (At — 0.1, original coordinates, 
supremum norm, Eq ~ 124 resulting in RR ~ 0.05). Upper bounds found obtained for vertices given in ascending order of 
the associated local phase space density (in comparison to the original (temporal) order of vertices) are H* = 2.71 x 10~ 4 
and 2.53 x 10 -4 for the Rossler and Lorenz system, respectively. (C,D) Dependence of the mean value and standard deviation 
(error bars) of the Hamming distance (obtained from 100 permutations) on the length TV of the underlying time series for (C) 
Rossler and (D) Lorenz system, suggesting that ~ TV -1 (dashed lines). 



For practical applications, Xu et al. p008 suggested that the possibility of classifying dynamical 



systems based on the motif distributions of adaptive nearest neighbor networks (see Sec. 4.1) is robust 
to variations in the choice of both Eq (the number of edges drawn from each node) and the motif order. 
It is certainly true that varying Eq (over a reasonable range) does not affect the corresponding results 
significantly. The choice of the motif order is, however, rather limited due to practical reasons — 2- and 
3-motifs offer very little scope, 5-motifs and higher orders are combinatorial nightmares and very quickly 
become computationally intractable. Hence, the choice of a metric based on 4-motifs is largely one of 
practical expedience. 
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3.2.2. 



-recurrence networks 



2010a|, whe re it has 
Gao fc Jin[ [2009a|b 



The problem of threshold selection has been discussed in detail by Donner et al. 
been shown that simple heuristics such as the turning point criterion proposed in 

(i.e., determining e by the - supposedly unique - turning point of the RR(e) relationship) can provide 
misleading results. Moreover, the thresholds proposed by such general (system-independent) criteria can 
depend crucially on both sampling (see Fig. [9]) and embedding. Although there is no universal threshold 
selection criterion, some general considerations help fixing e at an appropriate value. On the one hand, if 
e is too small, there are almost no recurrence points. Hence, the information contained in the e-recurrence 
network is rather limited. On the other hand, if e takes too large values (which is typically the case for 
the turning point criterion), every vertex is connected with many other vertices irrespective of their actual 
mutual proximity in phase space. One reasonable trade-off between these two extreme cases is choosing an e 
within the scaling region of the correlation integral, which coincides with the classical strategy for estimating 



the correlation dimension D2 using the Grassberger-Procaccia algorithm Grassberger & Procaccia 1983 



Following independent arguments, Schinkel et al. 20081 suggested selecting e for applications of RQA 



corresponding to recurrence rates RR < 0.05. Given sufficiently large N, this choice also provides reasonable 



local information about the attractor topology in phase space based on e- recurrence networks |Marwan 
Donner et al] |2010b|a] 



et al. , 2009 
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Fig. 9. Sampling effect on the dependence of the link density p of an £-recurrence network (equivalent to the recurrence rate 
RR) on e for (A) Lorenz system and (B) Rossler system with a = 0.2, b = 0.2 and c = 5.7 (N = 1,000, Euclidean norm, 
original coordinates). The three curves correspond to the sampling rates Ati = 0.01, At2 = 0.1, and At% — 1.0. Small circles 
indicate the recurrence thresholds suggested by the turning point criterion [Gao &: Jin| |2009a|b| . 



The quantitative characteristics of e-recurrence networks depend on e, which is of particular importance 
for global network measures. Specifically, the average path length C is approximately inversely proportional 
to e |Donner et al. 2010b| . For other measures such as transitivity T, global clustering coefficient C or 
assortativity 1Z (i.e., the correlation coefficient between the degrees of all pairs of directly connected 
vertices), the behavior varies with the system under study. As a general observation, for large e, C — > 1 
due to the increasing coverage of the attractor by the ^-neighborhoods. An approximate analytical theory 
for one-dimensional maps has been given by Donner et al. 2010b| . Corresponding statements hold for 



the transitivity T as well. The assortativity, however, shows a more diverse behavior: For e being small 
compared to the attractor diameter, we find a tendency towards smaller values as e increases, whereas for 
large £,7^—7-1 since k v — > N — 1 for all vertices. Similar observations can be made for vertex and edge 
properties. However, their spatial distributions usually remain qualitatively robust as long as e does not 
become too large |Donner et al. 2010a . For sufficiently large N, the features revealed by measures, such 
as centralities or local clustering coefficient, can be related to finer structures in phase space for small e, 
whereas there is a successive smoothing as the recurrence threshold increases. 

When comparing different time series from the same system, it is often desirable to fix the recurrence 
rate RR instead of e. Firstly, the resulting e-recurrence networks have approximately the same number of 
edges, which allows comparing the resulting topological properties of different networks more objectively. 
Secondly, the attractor diameter in phase space can change with varying control parameters. The question 
whether to apply a fixed RR or a fixed e is especially important when cases with rather different dynamical 
properties are to be compared (e.g., periodic and chaotic orbits), where the respective RR(e) relationships 
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are hardly comparable Zou et al. , subm. . 



4. Applications 



4.1 



Classification of dynamical systems 

Xu e£ a/. |2008| showed that the motif prevalence in adaptive nearest neighbor networks — in particular, 
the motif superfamily membership (i.e., the qualitative coincidence of the motif distributions within a 
large class of complex networks) — can be used to classify dynamics as chaos (with one positive Lyapunov 
exponent), hyperchaos (multiple positive Lyapunov exponents), noise, or a periodic orbit. As a real- world 



example, we apply the same method to experimental data (partially depicted in Fig. 10) of sustained tones 
voiced on a standard Bb clarinet over the dynamic range of the instrument — from E3 to Hq in standard 
scientific notation. The 20 distinct notes where individually recorded and manually preprocessed to extract a 
stationary (in terms of amplitude) period of data which was then smoothed and down-sampled. Specifically, 
each signal was recorded at 44.1kHz and consisted of 70,000 samples from the stationary sustained phase 
of the intonation. This was then down-sampled to a level with approximately 25 samples per cycle. The 
time delay was chosen to be the first minimum of mutual information (typically between 3 and 6) and the 
embedding dimension was 10. From the embedded time series, adaptive nearest neighbor networks have 
been constructed using the Euclidean norm and four neighbors of every vertex (Eq = 4). 




" n 

LL U 



0.005 0.01 0.015 
time (sec.) 



0.02 




0.005 0.01 0.015 
time (sec.) 



0.02 




0.005 0.01 0.015 
time (sec.) 



0.02 




0.005 0.01 0.015 
time (sec.) 



0.02 



10 




2000 4000 6000 8000 10000 
frequency (Hz) 



10 u 




2000 4000 6000 8000 10000 
frequency (Hz) 



10" 




2000 4000 6000 8000 10000 
frequency (Hz) 



10" 




2000 4000 6000 8000 10000 
frequency (Hz) 



Fig. 10. Bb clarinet tones, from E3 (below middle C) to Bg (above the treble clef), are analyzed with the adaptive nearest 
neighbor network method. In this figure we illustrate short sections of just three of these signals. On the left are short (20 
ms) samples of the sound wave for Bg, C5, F3 and E3 in the time domain. On the right we depict a corresponding sample 
power spectra for each. Note that the notes E3 and F3 are near the bottom of the lowest (chalumeau) register of the clarinet 
and have a characteristic rich and woody tone. The notes C5 and Bg are in the intermediate clarion register (where notes 
are produced by overblowing — which removes the lowest frequencies) and have a more pure, bright and penetrating sound. 
Notes in the altissimo register were not studied. 



Four representative networks are depicted in Fig. 11 Despite the striking variation in the appearance 
of these networks (each of which is characteristic of that particular recording) we find that the mesoscopic 
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network structures are remarkably similar. In particular, following [Xu et af. |2008 we compute the fre- 



quency of motif patterns for motifs of order 4. We find that 17 of the 20 distinct tones generate the same 
motif prevalence — these tones all belong to the same motif superfamily (see top line in Fig. 12). The 
remaining three tones belong to a distinct, but closely related family (bottom line in Fig. 12). 




Fig. 11. Adaptive nearest neighbor networks for four distinct tones on the clarinet. The lower two plots E3 and F3 correspond 
to the bottom of the Bb clarinet's range (in the top half of the bass clef), the note C5 is at the bottom of the clarion register 
(center of the treble clef) and Bg is at the top. 
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Fig. 12. Motif prevalence for the 20 analyzed clarinet tones (ordered according to their total frequency). Of these 20 tones, 
17 belong to the same motif superfamily — depicted in the upper row of this figure. The remaining three tones E3, F4 and 
Bg belong to a distinct, but closely related motif family — depicted in the lower row. 



While there is nothing concrete to link the three odd ball tones (E3, F4 and Bq), it is possible (due to 
the skill, or lack thereof of the musician — or the quality of the instrument]^]) that these three tones may 
be prone to less stationarity. The tone E3 is the lowest that the Bb clarinet will produce and the sound 
typically has more of a vibratory quality. The notes F4 and are in distinct registers of the instrument, 



A poor workman will blame his tools, in this case the intonation was performed by a poor musician on a cheap instrument. 
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but are produced with a similar length resonance (similar fingerings) either with or without overblowing. 
Hence, the dynamics of these tones should be mechanistically similar. Of course, they should also be similar 
to many other related tones. Nonetheless, the two motif superfamilies are very similar to one another (there 
is only one transposition). 

The motif super- family for these three mentioned notes (E3, F4 and Bg) and the remaining seventeen 
notes are distinct from all those reported in |Xu et a~. , 2008| . Nonetheless, both motif superfamilies are 
most similar to that of chaotic or hyperchaotic dynamics (each is only one permutation from the motif 
superfamily observed for hyperchaos). In all cases, the dynamics we observe in this experimental data is 
clearly distinct from that observed for a noisy periodic orbit. 

The indication of this application — that pure sustained clarinet tones are characteristically aperiodic, 
and consistent with chaotic or hyperchaotic dynamics — is intriguing, but also preliminary. Further work 
with both network based methods and other techniques from nonlinear time series analysis is required. 



4.2. Identification of dynamical transitions 

One of the major applications of traditional RQA is the identification of dynamical transitions from time 
series. The RQA measures can be calculated for small square windows of size w moving along the main 



k+w- 
*J I i,j=k 



Trulla et al. 



1996 



Marwan et al. 



2007b I . This approach 



diagonal of the RP, i.e., in the sub-RP R 

allows studying the temporal variation of the different KQA measures, and, hence, identifying transitions 
in the dynamics of the studied system in terms of significant changes of these measures with time. For 
example, it has been shown that the diagonal line-based RQA measures are able to detect transitions 
between chaotic and regular dynamics in maps, whereas vertical line-based measures can identify chaos- 
chaos transitions 



phases. 



Marwan et al., 2002b , e.g., in terms of detecting different properties of the laminar 



wan et al. 



Similar results have been obtained for the quantitative characteristics of e-recurrence networks. Mar- 
~ |2009| studied the bifurcat ion cascade of the logistic map x n+ \ = ax n (l — x n ) using both 
RQA measures and global properties of the recurrence networks associated with individual realizations 
for different values of a. It has been found that the presence of periodic windows is clearly detected by 
both transitivity (and global clustering coefficient) as well as average path length. Specifically, periodic 
dynamics is indicated by T = 1 (C = 1) and C = 1, whereas, for example, the recurrence rate still shows 
different values in dependence on the specific period. In addition, it has been observed that for finite 
e, sudden jumps of C precede band merging points due to a merging of formerly disconnected network 
clusters. Additional pronounced minima of C have been found to coincide with chaos-chaos transitions. 
Although these bifurcations are also detectable with vertical line-based RQA measures such as laminarity 
or trapping time |Marwan et al. , 2002b| , the shifts in the average path length are particularly well localized 
at the appropriate values of a. A similar study of complex bifurcations in a two-dimensional parameter 
space of the time-continuous Rossler system has recently been reported by Zou et al. subm.| . A comparison 
with maximum Lyapunov exponents obtained from long realizations of the system revealed that recurrence 
network measures estimated from short time series allow a reasonable distinction between periodic and 
chaotic windows, which performs somewhat better than a corresponding discrimination based on RQA 
measures. 

As a real- world example for the detection of hidden transitions by means of e-recurrence networks, we 
reconsider the analysis of a marine terrigenous dust flux record from the Ocean Drilling Program (ODP) site 
659 (Tiedemann et aZ. , |1994] (see Fig. 134), which is located in the Atlantic close to Northwest Africa. This 
time series has a length of iV = 1,221, covering the last 5.0 Ma (million years) with an average sampling 
time of 4.1 ka (thousand years). Note that the time scale is not equidistant, the standard deviation of 
sampling time being 2.7 ka. Compared to the long geological time span covered, this deviation can be 
considered as quite small. 

The ODP 659 terrigenous dust flux record has been used to infer epochs of arid continental climate 
conditions and related long-term changes in the African climate. Various studies based on this record 
have restricted themselves to the use of linear methods of time series analysis [deMenocal 1995; Ravelo 



et al. , 2004 ; Trauth et al. , 2009 , revealing changes in dominating cyclic (Milankovich) components and 
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Fig. 13. (A) Terrigenous dust flux record of ODP site 659, and corresponding network measures (B) average path length C, 
(C) transitivity T, (D) assortativity coefficient 1Z and (E) diameter D obtained from e-recurrence networks (m = 3, r = 2, 
p = 0.05, supremum norm, window size of 410 ka with 90% overlap). The dotted vertical lines indicate time intervals that are 
identified as marked features with respect to the simple confidence intervals described in the text, the dash-dotted horizontal 
lines correspond to the mean values for the null-model. 



their possible relationship with known globally observable climate shifts such as the onset of Northern 
hemisphere glaciation, the mid-Pleistocene climate shift, or the intensification of the Walker circulation 
after about one million years before present (BP) |Mudelsee fe Raymo 2005 St. John & Krissek 2002 



McClymont & Rosell-Mele 


2005 . 


Recently, 


Marwan et al. 


|2009 



studied this time series by means of e-recurrence networks. Note again 
that unlike most other existing methods of time series analysis, RP based techniques do not explicitly 
require a regular sampling. The only implicit assumption is that the data used in this approach repre- 
sents the distribution of observations in the underlying phase space in a statistically reasonable way. This 
presumption makes RPs and recurrence networks natural candidates for the investigation of paleoclimate 
time series, since irregular sampling is a typical problem in the analysis of this kind of data. 

In the following, we will reconsider the recent findings and provide complementary results from ad- 
ditional network characteristics, using a time series that extends about 500 ka further back into the past 



than that used by Marwan et al. 



2009] . For consistency, we apply time-delay embedding with the optimum 
2 (selected based on the false nearest neighbor and mutual informa- 

0.05 



embedding parameters m = 3 and t 

tion criteria) and a variable recurrence threshold e conserving a constant RR 



Schinkel et al. 2008 
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Donner et al., 2010a| (cf. Sec. 3.2.2). To study transitions in the dust record, we construct e-recurrence 



networks for 112 windows of size 100 samples (corresponding to periods of on average 410 ka) covering the 
line of identity (i = j) with a mutual overlap of 90% (resulting in a step size of 10 samples, or approx. 
41 ka). To determine the time scale of the windowed network measures, we chose the windows' mid-points. 

We furthermore statistically test whether the network characteristics at a certain time differ signifi- 
cantly from the general network characteristics expected given the phase space distribution of state vectors 
for the whole embedded marine dust record and chosen window size. The corresponding null- hypothesis 
is that the network measures observed for a certain window are consistent with being calculated from a 
random draw of 100 state vectors from the prescribed phase space distribution induced by the entire time 
series. We can justly assume a thus randomized embedded time series without loosing essential informa- 
tion, because network measures are permutation-invariant (a similar test for RQA measures requires a 
more advanced method [Schinkel et al. 2009| ). In order to create an appropriate null- model, we use the 
following approach: 

(i) Randomly select w = 100 state vectors x a from the complete embedded time series. Here, the specific 
choice of w corresponds to the chosen window size. 

(ii) Use this random sample of state vectors for constructing an e-recurrence network. 

(iii) Calculate the network measures of interest from this recurrence network. 

(iv) Repeating this procedure 50, 000 times, we obtain a test distribution for each of the network measures. 
The 5% and 95% quantiles of the true test distribution, which can be estimated from these empirical 
distributions with sufficiently high confidence, can then be interpreted as 90% confidence bounds (see 
Fig.pfe-E). 



With this approach, one may test whether the spatial distribution of state vectors obtained for a given time 
slice is typical for the whole time series. Hence, the suggested procedure tests in fact against stationarity of 
certain geometric phase-space properties of the system under study. Time intervals yielding values of some 
network property that significantly differ from the corresponding distribution obtained from the random 
samples can be interpreted as possibly containing changes in the phase space structure and, hence, the 
observed dynamics encoded in the considered time series. 

The network measures average shortest path length C and transitivity T exhibit a distinct variability 
(Fig. [13^3 and C). As the most remarkable features, C highlights epochs of significantly increased values 
between 3.5 and 3.3 Ma, around 2.1, and between 1.9 and 1.7 Ma BF(^} T discloses epochs of increased 
values between 3.5 and 3.0 Ma as well as between 2.5 and 2.0 Ma. The assortativity coefficient 1Z as a 
measure of the continuity of the density of states |Donner et al. , 2010b| is expected to show some correlation 



with transitivity and global clustering coefficient. In fact, 7Z also increases significantly between 3.5 and 
3.2 Ma BP (Fig. [l3p). The evolution of 1Z is, however, more distinct from T after ~3.0 Ma BP, where 7Z 
decreases markedly especially between 0.9 and 0.5 Ma BP. Moreover, there appears to be a slight trend in 
the evolution of 1Z, resulting in a tendency to decrease from higher values in the distant past towards lower 
values in the present. The network diameter D (i.e., the maximum shortest graph distance between all 
pairs of vertices) evolves similarly to the average path length C (Fig. [13^), since both measures quantify 
statistical properties of the distribution of shortest path lengths on the network. Note that the relative 
amplitudes of D during the epochs of significant increase differ markedly from those of C. 

The time intervals identified by the different complex network measures are robust and seem to be well 
correlated with some major transitions in the climate system (e.g., the end of the Pliocene optimum at 



about 3.4-3.1 Ma BP) [Marwan et al. 2009]. We note that these specific intervals have not yet been found 
using classical methods of time series analysis, such as spectral analysis or breakpoint regression Trauthl 
et al. , 2009 , and do also slightly differ from recent results obtained using RQA |Marwan et al. 



We relate this to the fact that recurrence network characteristics indeed capture conceptionally different 



2008 



structural properties of a dynamical system from rather short time series (see Sec. 2.3.4). A detailed 
climatological interpretation of our findings will be given in an upcoming paper. 



4 In paleoclimatology, BP (before present) refers to the number of years before the reference year 1950. 
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5. Summary 

Recurrence is a fundamental property of many dynamical processes. It is a concept successfully used 
in the study of dynamical and complex systems, and time series analysis, e.g., using recurrence times 
statistics, first return maps, recurrence plots, or recurrence quantification analysis. Conversely, network 
theory provides important insights in the study of many complex systems. By exploiting the duality 
between the recurrence matrix in the study of dynamical systems and the adjacency matrix of a complex 
network, we have demonstrated how information about dynamical recurrences can be used to construct 
complex networks from time series. These recurrence-based complex networks provide a new approach for 
time series analysis and offer a promising and complementary view for the study of dynamical systems. 
Applying well established complex network measures, we are able to characterize and classify the dynamics 
of complex systems, to detect dynamical transitions, or identify invariant substructures. 

The quantitative characteristics of recurrence-based complex networks have a clear interpretation in 



terms of geometric properties of the underlying system in phase space (Donner et al., 2010b] . In partic 



ular, recurrence networks do not take time information into account and, hence, do not explicitly rely 
on the presence of equally spaced observations, which is an important problem for the analysis of many 
real- world time series. The only implicit assumption one has to make is that the actual phase space density 
of the system is sufficiently well represented by the given set of points. Specifically, temporal correlations 
between individual observations are not taken into account, which makes recurrence-based complex net- 
works distinctively different from the majority of other methods of time series analysis. Sufficiently long 
realizations guarantee stable network structures that do not depend on the specific realization of the sys- 
tem. However, besides the loss of information about temporal structures, the purely geometric point of 
view on higher-order statistical properties of a system is partially related with higher computational costs 
for estimating certain complex network measures {e.g., betweenness centrality) in comparison with other 
methods. Furthermore, it should be noted that there are alternative methods of complex network-based 
time series analysis (e.g., transition networks, visibility graphs, or correlation networks, cf. Sec. [2]), which 
are not based on the recurrence concept and, hence, are characterized by distinctively different properties. 

The new approach of recurrence-based complex networks combines two successful concepts in modern 
complex systems studies: the recurrence plot and the complex network. The first promising applications 
of this approach illustrate the potential of recurrence networks and their interdisciplinary relevance. The 
underlying conceptual idea of constructing networks on the basis of mutual proximity relations in phase 
space is rather simple and thus has (as RQA Webber Jr. et al, 2009] ) the potential to be applied in a 



meaningful way in various fields of science. We have to conclude, however, that the new method of recurrence 
network analysis is still in its infancy. A specific question that is necessary to be systematically addressed is 
which specific approach to prefer for which particular application. For example, all previous applications to 
well-known paradigmatic models considered only relatively low- dimensional systems, whereas recurrence 
network properties have not yet been explicitly studied for high-dimensional systems. We emphasize that 
the problem of appropriate embedding and the available amount of data can be expected to become 
more crucial as the dimension of the system under study increases. In summary, there are many open 
questions concerning the specific features and applicability of this new conceptual approach, which will 
surely stimulate future investigations. This statement also holds for the other recent approaches to analyzing 
time series by means of complex network methods. 
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